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We calculate the properties of neutron matter and highlight the physics of chiral three-nucleon 
forces. For neutrons, only the long-range 2-7r-exchange interactions of the leading chiral three- 
nucleon forces contribute, and we derive density-dependent two-body interactions by summing the 
third particle over occupied states in the Fermi sea. Our results for the energy suggest that neutron 
matter is perturbative at nuclear densities. We study in detail the theoretical uncertainties of the 
neutron matter energy, provide constraints for the symmetry energy and its density dependence, 
and explore the impact of chiral three-nucleon forces on the S-wave superfluid pairing gap. 
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I. INTRODUCTION 

The physics of neutron matter ranges over excit- 
ing extremes: from universal properties at low densi- 
ties P, Q that can be probed in experiments with ul- 
tracold atoms to using neutron matter properties at 
nuclear densities to guide the development of a universal 
density functional 0, H| and to constrain the physics of 
neutron-rich nuclei; to higher densities involved in the 
structure of neutron stars [bL In the theory of nuclear 
matter, recent advances 07HI are based on systematic 
chiral effective field theory (EFT) interactions [t| [l(| 
combined with a renormahzation group (RG) evolution 
to low momenta [HI, (Hf. This evolution improves the 
convergence of many-body calculations 0, E3, EU and the 
nuclear matter energy shows saturation with controlled 
uncertainties [1[. In this paper, we extend theses devel- 
opments to neutron matter with a focus on three-nucleon 
(3N) forces. 

Our studies are based on evolved nucleon-nucleon 
(NN) interactions at next-to-next-to-next-to-leading or- 
der (N 3 LO) [H EH and on the leading N 2 LO 3N 
forces [13, EH- I n Sect. [TlJ we show that only the long- 
range 27r-exchange 3N interactions contribute in pure 
neutron matter. We then construct density-dependent 
two-body interactions V 3 n by summing the third par- 
ticle over occupied states in the Fermi sea. Effective 
interactions of this sort have been studied in the past 
using 3N potential models and approximate treatments 
(see for example, Refs. [IiJ[2(|). We derive a general op- 
erator and momentum structure of V3N and analyze the 
partial- wave contributions and the density dependence of 
^3N- This provides insights to the role of chiral 3N forces 
in neutron matter. 

In Sect. IIII) we apply V3N to calculate the properties 
of neutron matter as a function of Fermi momentum k-p 
(or the density p = kp/(3n 2 )) based on a loop expan- 
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sion around the Hartree-Fock energy. Our second-order 
results for the energy suggest that neutron matter is per- 
turbative at nuclear densities, where N 2 LO 3N forces pro- 
vide a repulsive contribution. We study in detail the the- 
oretical uncertainties of the neutron matter energy and 
find that the uncertainty in the C3 coefficient of 3N forces 
dominates. Other recent neutron matter calculations lie 
within the resulting energy band. In addition, the energy 
band provides constraints for the symmetry energy and 
its density dependence. Finally, we study the impact of 
chiral 3N forces on the superfluid pairing gap at the 
BCS level. We conclude and give an outlook in Sect. [TV] 



II. THREE-NUCLEON FORCES AS 
DENSITY-DEPENDENT TWO-BODY 
INTERACTIONS 

Nuclear forces depend on a resolution scale, which is 
generally determined by a momentum cutoff A, and are 
given by an effective theory for scale-dependent two- and 
corresponding many-nucleon interactions j9Mlll |21| : 

H(A) = T + Vnn(A) + V3N (A) + V4n(A) + .. . . (1) 

Our calculations are based on chiral EFT interactions. 
We start from the N 3 LO NN potential (A = 500 Me V) 
of Ref. [IH and use the RG to evolve this NN poten- 
tial to low- momentum interactions Viowfc with a smooth 
n cxp = 4 regulator with A = 1.8 - 2.8 fm" 1 [HH^. This 
evolution softens the short-range repulsion and short- 
range tensor components of the initial chiral interac- 
tion 0, [23| . Based on the universality of Vy ow k [1, EI| > 
we do not expect large differences starting from different 
N 3 LO potentials. 

In chiral EFT without explicit Deltas, 3N forces start 
at N 2 LO and contain a long-range 27r-exchange part V c , 
an intermediate-range l7r-exchange part Vd and a short- 
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The 27r-exchange interaction is given by 
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where = — denotes the difference of initial and 
final nucleon momenta (i,j and k = 1,2, 3) and 
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while the l7r-exchange and contact interactions are given 
respectively by 
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Ve = zT^r E( t j • rfc ) 



( Ti • Tj) (<Ti ■ qj) , (5) 
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with g A = 1.29, / w = 92.4 MeV, = 138.04 MeV and 
A x = 700 MeV. For 3N interactions, we use a smooth 
regulator as in Ref. [8j, 

(p 2 + 3</ 2 /4) 2 



fn(p,q) = exp 
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with a 3N cutoff A 3 nf that is allowed to vary indepen- 
dently of the NN cutoff and probes short-range three- 
body physics. Here, p and q are initial Jacobi momenta 
(the final Jacobi momenta are denoted by p' and q'). The 
exchange terms of the 3N force are included by means of 
the antisymmetrizer 

A 123 = (1 + P12P23 + ^i3-P 23 )(l - P23) , (8) 

= 1 — P\1 — P\3 — P-23 + P12P23 + P13P23 7 (9) 

where Pjj is the exchange operator for spin, isospin and 
momenta of nucleons i and j. The regulator /r(p, is to- 
tally symmetric when expressed in the nucleon momenta 
ki, and thus the direct and exchange terms contain the 
same regulator. 

In this work, we take the N 2 LO 3N forces as a 
truncated basis for low-momentum 3N interactions and 
assume that the c$ coefficients of the long-range 2tt- 
cxchange part V c are not modified by the RG evolu- 
tion. This follows the strategy adapted in Refs. 0,11, H3], 
until it will be possible to evolve many-body forces in 



momentum space starting from chiral EFT. For chiral 
low-momentum interactions, the cd and ce couplings 
have been fit for various cutoffs to the 3 H binding en- 
ergy and the 4 He radius in Ref. Q. However, as will be 
shown in the following, in pure neutron matter only the 
Ci and c 3 terms of the 27r-exchange part V c contribute, 
so that the leading low-momentum three-neutron inter- 
actions will not depend on the shorter-range parts. Since 
the cd and ce couplings do not enter in neutron mat- 
ter, we take for the Cj coefficients the consistent values 
used in the N 3 LO NN potential of Ref. [HI , in particular 
ci = -0.81 GeV -1 and c 3 = -3.2 GeV -1 . 

In addition, we will estimate the uncertainties of the 
two assumptions we have made for 3N forces. First, we 
will vary the cutoff to probe the sensitivity to neglected 
short-range many-body interactions and to the complete- 
ness of the many-body calculation. Second, we will study 
the dependence of our results on the choice of the c, coeffi- 
cients within their theoretical uncertainties in Sect. MI "Dl 



A. Chiral 3N forces in neutron matter 

Neutron matter presents a very interesting system, be- 
cause only certain parts of the N 2 LO 3N forces contribute 
(see, for instance, Ref. 25]). First, the contact inter- 
action Ve vanishes between antisymmetrized states, be- 
cause three neutrons cannot interact in relative S-wave 
states due to the Pauli principle. 

Second, the antisymmetrized l7r-exchange part Vd 
vanishes due to the particular spin-momentum structure 
of this interaction. With r, • t< = 1 for neutrons, one has 
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where in Eq. (|10[) and in the following the isospin ex- 
change operator PJj = 1 in the antisymmetrizer. The 
sum over i =/= j =/= k can be grouped into three terms, 
<ti • qi (er 2 + er 3 ) • q 1; and with 1 (23) replaced by 2 (13) 
and 3 (12). The first of these three terms is independent 
of the momenta of particles 2 and 3 and is nonvanishing 
only if the spin part of the wave function is symmetric 
under the exchange of the spins of particles 2 and 3. This 
implies P 2 3 = 1 and using the product representation of 
the antisymmetrizer, Eq. ([5]). leads to 



^123 o"i ■ Qi (o"2 + 03) • Qi = . 
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The other two terms in Vd vanish similarly. The physical 
reason for A123VD | = is that the two particles inter- 
acting through the contact interaction cd of Eq. ([2]) are 
required to be in a symmetric two-body spin state due 
the structure <Xj + (Tj (because pion exchange couples to 
the sum of the spins), but the interaction does not de- 
pend on the momenta of these two particles, which means 
the momentum-space wave function will be symmetric. 
Therefore, the two particles cannot be in an overall anti- 
symmetric state. 
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FIG. 1: Diagrammatic representation of the antisymmetrized N 2 LO 3N force in neutron matter. The C4 term in V c and the 
shorter-range parts Vd and Vb vanish in neutron matter. 



Finally, the C4 term in Eq. (Q} of the 27r-exchange 
part V c does not contribute for neutrons, because of the 
isospin structure with (nnn\Ti ■ (T2 x T^)\nnn) — 0. As a 
result, the N 2 LO 3N force in neutron matter is given by 
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The direct, single- and double-exchange terms included 
in Eq. (fj"2"]) are shown diagrammatically in Fig. [T] Since 
C3 is typically a factor 4 larger than c\ and because qi • qj 
increases with density compared to m 2 (for densities of 
interest q ~ fcp ~ 2m 7r ), the N 2 LO 3N force in neutron 
matter is therefore dominated by the C3 contribution. 



B. Operator structure of V3N 

Based on the N 2 LO 3N force, Eq. (fl"2")) . we construct 
an antisymmetrized density-dependent two-body interac- 
tion 1^3N in neutron matter by summing the third particle 
over occupied states in the Fermi sea, 
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where the trace is over the spin of the third particle and 
Tit denotes the Fermi-Dirac distribution function at zero 
temperature, rik = 9(kp — |k|). The antisymmetrization 
in Eq. (|13p can also be written in a more symmetric form 
using normalized and antisymmetrized three-body states, 
|123) as = ^^li 23 1 123), so that (T2'3'| A 23 ^3n|123) = 

as (l'2'3'|V A 3 N |123) as . Diagrammatically, Eq. (33]) corre- 
sponds to connecting one incoming and one outgoing line 
in the antisymmetrized 3N force with a noninteracting 
propagator: 
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The resulting density-dependent two-body interaction 
V3N corresponds to the normal-ordered two-body part of 
3N forces (for normal- ordering, see Refs. (2(| H3|), where 



the normal ordering is with respect to the Fermi sea. As 
a result, the antisymmetrized V3N is added to the anti- 
symmetrized NN interaction Vnn, as = (1— P12) Vnn when 
evaluating two-body contributions beyond the Hartree- 
Fock level. For the Hartree-Fock (or normal-ordered 
zero-body) contribution, the 3N diagram has a symme- 
try factor 1/6, so that I^3n/3 is added to Vnn, as at the 
Hartree-Fock level (see Sect. IIII A[) . For one-body con- 
tributions, normal ordering leads to V^3n/2 being added 
to Vnn, as, for example, for the first-order contribution 
to the self-energy. These factors ensure the correct sym- 
metry factors for diagrams involving NN and 3N interac- 
tions. The above normal-ordering factors may not always 
be taken into account properly in the literature [2js-l30| 
(for example, it is incorrect to use Vnn, as + V^n/S or 
VNN,as+^3N for both Hartree-Fock and higher-order con- 
tributions). 

In addition to the density dependence, the two-body 
interaction V3N depends on the incoming and outgoing 
relative momenta, k = (ki — k2)/2 and k' = (k^ — k 2 )/2, 
on the spin of particles 1 and 2, and on the two-body 
center-of-mass momentum P = ki + k 2 = + k 2 . How- 
ever, as will be shown in Sect. IIII A[ the P dependence 
is weak compared to the dependence on the relative mo- 
menta. This can be understood from the momentum 
transfers in the pion propagators of Fig. [TJ The momen- 
tum transfers in the first two diagrams are k — k', (for 
the left and right propagators) and k + k', 0. While these 
terms vanish in the considered case of V c , the only weak 
P dependence would arise from the regulator functions 
in Eq. (|12p. For the remaining four diagrams in Fig. [T] 
the momentum transfers are k ± k' , P/2 ± k' — k3 and 
P/2+k-k 3 , P/2±k'-k 3 . Therefore, the P dependence 
is only through kg = P/2 — k3 and can be transformed 
into a dependence due to rip /2-k' (by substitution in the 
integral) and through the regulator functions. For typi- 
cal center-of-mass momenta P < fcp, this results in the 
weak P dependence. 

For the derivation of V3N, we use Mathematica and 
automate the momentum and spin exchange operations 
by representing all spin operators in matrix form. The 
antisymmetrized 3N force in neutron matter, Eq. (|12[) . is 
represented in this explicit three-particle spin basis and 
for general particle momenta k^. We then perform the 
spin trace over the third particle and project the resulting 
two-body matrix on a complete set of two-body spin op- 
erators. In this way, we obtain a general operator struc- 
ture of V^3N in spin-saturated neutron matter in terms of 
a set of basic integral functions. Only at this point we 
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transform V^n to the relative and center-of-mass basis 
and adapt a fixed- P approximation. Since the number of 
spin operators is minimal for P — 0, we take P = for 
simplicity and our final result is given by 



where the functions Ai 2 (k, k') and Pi 2 (k, k') include all 
spin dependences, 
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A 12 (k, k') = 2 p 2 (k, k') + 2a\k, k') - a\k, -k') - 6 (k, k') 
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S 12 (k + k') p° (k, k') - S 12 (k - k') p°_(k, k') -1*1*1 d ab (k, k') - d ab (k, -k' 



2i(<r 1 + < r 2 )» e°(k,k')-c°(k,-k') 



P 12 (k, k') = -2 p* (k, k') + 2a 2 (k, k') - a 2 (k, -k') - b (k, k') 
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and the basic integral functions are defined by 
(k ± k')" 



d ab (k,k')-d ab (k,-k') 



p n ± (k } k') 
a n (k,k') 
b n (k,k') 
c"(k,k') 

C & (k,k') 



((k±k') 2 



1 



((k + k 3 )-(k'+k 3 )) n 



k3 ((k + k 3 )2 + m2)((k' + k 3 )2 + m2) ' 

((k + kQ.(k + k 3 ))" 
k3 ((k + k') 2 + m2)((k + k 3 )2+m2) ' 

((k + k 3 ) ■ (k' + k 3 ))"((k + k 3 ) x (k' + k 3 )) a 

((k + k 3 ) 2 + ml) ((k' + k 3 ) 2 + ml) 
(f1r ,^ n, , N ,» (k + kQ a (k + k 3 ) b + (k + kO b (k + k 3 ) a -|,5 ab (k + kO.(k + k 3 ) 
k3 U j ' 1 3jj 2 ((k + k') 2 + m 2 ) ((k + k 3 ) 2 + m 2 ) 



r 



(16) 



(17) 

(18) 
(19) 
(20) 
(21) 
(22) 



In Eqs. (|16|) (|22|) . the indices a, b run over the three com- 
ponents of the spin operators, the tensor operator is given 
by Si 2 (p) = (cri ■ p)(*2 ■ p) -p 2 *i ■ c 2 /3, and the over- 
line denotes a symmetrization in the relative momentum 
variables, x(k,k') = x(k, k') -I- x(k',k). In addition, we 
have introduced the short-hand notation for the integra- 
tion over the momentum of the third particle, 



dka 



n k3 f R (k',k 3 )f R (k,k 3 ), (23) 



/ k3 J (2^) 3 

where the regulator, Eq. ([7]), expressed in terms of the 



relative and third-particle momenta for P — is given 
by f R (k, k 3 ) = exp[-(fc 2 + fc 2 /3) 2 /A4 NF ]. 

The density-dependent two-body interaction V 3 n in- 
cludes all spin structures that are invariant under com- 
bined rotations in spin and space in a spin-saturated sys- 
tem. In addition to the central spin-independent and 
spin-spin (cri ■ er 2 ) interactions, the functions Ai 2 and 
B\ 2 include tensor forces (Si 2 ), spin-orbit interactions 
(c™-terms) and additional tensor structures terms), 
which can be expressed in terms of S\ 2 and quadratic 
spin-orbit interactions. Finally, we have checked that 
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FIG. 2: (Color online) Diagonal momentum-space matrix ele- 
ments of the density-dependent two-body interaction V3N for 
P = in the 1 So channel. Results with A3NP = 2.0 fm -1 are 
shown versus relative momentum k for different Fermi mo- 
menta fcp = 1.0, 1.2, 1.4 and 1.6 fm - (increasing in strength) 
in neutron matter. 



^3N is even (odd) in the scattering angle 0k, k' in the 
two-body spin S = (S = 1) channel. 



C. Partial- wave matrix elements of Vsn 

Next, we expand T^3N in two-body partial waves and 
show the diagonal momentum-space matrix elements in 
Figs. [2] and [3] in the 1 So and the spin-triplet P-wave 
channels for different Fermi momenta kp = 1.0,1.2,1.4 
and 1.6 fm - (in units where h 2 /m — 1, with nucleon 
mass m). The partial- wave matrix elements are nor- 
malized to the direct term (as for NN interactions). In 
Figs. [2] and [3l we have interpolated between the mo- 
mentum mesh points. To set the scale, a typical S- 
wave strength for low-momentum NN interactions is 
Mowfc(0,0) ~ -2.0 fm [O, 03, but the relative 3N 
strength depends on whether V^n is used at the Hartree- 
Fock, one- or two-body level (see Sects. MI Al to IlII C|) . 
We find that the dominant partial wave is the chan- 
nel, which is found to be repulsive. In addition, for 
kp = 1.6 fm - we plot in Fig. [3] the central and central 
plus tensor parts of V3N • Because Si 2 (p)/p 2 = —4/3,2/3 
and —2/15 in the 3 Po, 3 Pi and 3 P 2 channels, we conclude 
from Fig. [3] that the tensor force in V3N is repulsive. 
Moreover, with L • S = -2, -1 and 1 in the 3 P , 3 Pi 
and 3 P 2 channels, we observe that the spin-orbit inter- 
action in ^3N (determined by the difference of the solid 
and dashed lines) is attractive. 

The density dependence of the two-body matrix ele- 
ments depends on the partial wave. In the chan- 
nel this dependence can be approximately parameterized 
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FIG. 3: (Color online) Diagonal momentum-space matrix el- 
ements of the density-dependent two-body interaction I/3N 
for P = in the spin-triplet P-wave channels. Results with 
A3NF = 2.0 fm -1 are shown versus relative momentum k for 
different Fermi momenta kF = 1.0,1.2,1.4 and 1.6 fm - (in- 
creasing in strength). For kF = 1.6 fin , the dotted lines 
represent the central parts (degenerate in J) of V3N, whereas 
the dashed lines include the central plus tensor interactions 
(without the c" -terms in Eqs. (fl"6|) and fT7)l). 



by a power of the Fermi momentum, V3jsi(k',k,kp) « 
F{k' , k) {kp /kp) , with some reference Fermi momentum 
kp . In the spin-triplet channels this is more complex due 
to the different momentum and density dependences of 
the various operator structures in V3N- 

Finally, it is possible to improve the P = approxi- 
mation and to perform the sum over the third particle 
self-consistently, which corresponds to closing the line in 
Eq. (|14[) with the self-consistent propagator. For finite P, 
one has more allowed spin operators with more complex 
integral functions that can depend on P but also on the 
angle of P with respect to k and k'. Since 1^3N has been 



derived using Mathematica for general particle momenta 
kj, this is directly possible. One could then explore angle- 
averaging over P or averaging over the magnitude of P. 
However, as will be shown in the next section, the P = 
approximation is reliable for bulk properties and neu- 
tron matter based on chiral low-momentum interactions 
is sufficiently perturbative, which justifies using the non- 
interacting density to sum over the third particle. 



III. RESULTS 

We apply the developed density-dependent two-body 
interaction V3N to calculate the properties of neutron 
matter in a loop expansion around the Hartree-Fock 
energy. These are the first results for neutron matter 
based on chiral EFT interactions and including N 2 LO 3N 
forces. The many-body calculation follows the strategy 
of Refs. 0, H, HH, but with significant improvements for 
the second-order contributions involving F3N and with 
fully self-consistent single-particle energies. 



A. Hartree-Fock and P dependence of V3N 
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The contributions to the Hartree-Fock energy arc 
shown diagrammatically in Fig. |4] and the first-order NN 
and 3N interaction energies are given by 
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where V is the volume and we use the shorthand nota- 
tion i = kj<Xj in the bra and ket states. The momentum- 
conserving delta functions are not included in the NN 
and 3N matrix elements. It is evident from Eq. (j2"5| 
that the correct 3N symmetry factor is obtained when 
the antisymmetrized two-body interaction = (1 — 
P12) Viowfc + ^3n/3 is used at the Hartree-Fock level. 
With the expansion in two-body partial waves, we have 
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In Eq. (|26l) and in the following, all partial-wave ma- 
trix elements are normalized to the direct term, so that 



FIG. 4: Top row: Diagrams contributing to the Hartree-Fock 
energy. These include the kinetic energy -Ekin and the first- 
order NN and 3N interaction energies, Egl and E$. Middle 
and bottom rows: Second-order contributions to the energy 
due to NN-NN interactions E[ 2) , NN-3N and 3N-3N interac- 
tions, where 3N forces enter as density-dependent two-body 

respectively, and the remaining 

3N-3N diagram K~ 



interactions, E^l and E± 
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For N 2 LO 3N forces, the first-order 3N interaction en- 
ergy has been calculated exactly in Ref. (2|| , without the 
fixed-P approximation in V3N, 
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with kjj = ki — kj. By comparing the Hartree-Fock 
energy with the density-dependent two-body interaction 
V~ 3N for P = 0, Eq. {26]), to the result with the full 
first-order 3N energy, Eq. ([27j) . we can assess the reli- 
ability of the P = approximation in V~3N- Our re- 



sults for £^ +3N ff 



= E, 



kin 



E 



(i) 

NN 



E^ m comparison 



to the full Hartree-Fock energy E 



(i) 

NN+3N,full 



are shown 



for densities p < 0.16 fm (below saturation density 
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0.05 0.10 0.15 

P [fm 3 ] 

FIG. 5: (Color online) Energy per particle E/N of neutron 
matter as a function of density p at the Hartree-Fock level, 
based on low-momentum chiral NN and 3N interactions with 
A/A3NF = 2.0 fm" 1 . The first-order energy E^ eB using V3N 
with P — is compared to the exact Hartree-Fock energy 
i?3N full . To show the impact of chiral 3N forces, we also give 
the first-order energy based only on NN interactions. 

po = 0.16 fm~ 3 ) in Fig. [5] The P = approximation 
for V3N leads to energies that agree very well with the 
exact result at the Hartree-Fock level. For saturation 
density (fep = 1.7 fm -1 ), the difference is 340 keV, while 
for lower fcp, the differences are always below 100 keV, 
as documented in Table [I] This is very encouraging for 
evaluating higher-order contributions, where calculations 
with full 3N forces become increasingly more complex 
than with two-body interactions. 

As expected from the partial-wave matrix elements, we 



find in Fig. [5] that chiral 3N forces at first order increase 
the neutron matter energy with increasing density, com- 
pared to results based only on low-momentum NN in- 
teractions. This repulsion is due to the central parts of 
1^3N, because noncentral forces do not contribute to the 
Hartree-Fock energy. 



B. Second order: single- particle energies 

For the contributions beyond Hartree-Fock, we need 
to calculate the single-particle energies £ k . These are 
determined by the self-consistent solution to the Dyson 
equation. To second order, we have 

e k = ^ + £«(fc)+Rc£ (2) (Mk). (28) 
2m 

Here denotes the first-order NN and 3N contribu- 
tions to the self-energy, which are real, S'- 2 ' includes 
the second-order two-particle-one-hole and two-hole- 
one-particle terms, and the self-energy is averaged over 
the single-particle spin. The diagrams contributing to 
the self-energy to second order are given by 




The resulting first- and second-order self-energies can be 
expressed in terms of two-body partial waves and one has 
(see also Ref. (3l| ) 



£ (1) (fci) = ^ f k 2 2 dk 2 f d cos 9u u u 2 n U2 J2 (2J + l)<fci2/2|^ ) J |fc 12 /2)(l-(-l)'+ s+1 ), (30) 

S,l,J 

EWfo.wi) = / PdP I kdk I dk ' E E E iV ~ l " Y *Jfo Mk) IW(k') Y* m „(k') 

1 S,M S ,M' l,l',l",m,m',m" J,M 



x (C{™SM S ) 2 C?,% SM , C?Z»sm> s (k\V^j\k') (k>\V$ u \k) (1 - (-l)'^ 1 ) (1 - (-1)''+^) 
(l-n2 +k ,)(l-rc.£_ k ,)np_ kl ne +k , ri£_ k , (1 - n P _ kl ) 



(31) 



Wl + Ep_ kl - £p +k , - £p_ k , +iS Wi + E P _ kl - £p +k , - £p_ k , - id 

I 

where CfJf SM denote Clebsch-Gordan coefficients and axis is taken in the P direction, we have used kdk = 
Yi m (k) are spherical harmonics. In Eq. ([31]), the z- p ki dki/(4?r), and k = ki - P/2 determines the ar- 
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FIG. 6: (Color online) Effective mass m*(fcp)/m at the Fermi 
surface as a function of Fermi momentum &f in neutron mat- 
ter. Results for A/A3NF = 2.0 fm -1 are shown at the Hartree- 
Fock level, plus second-order contributions, and based only on 
NN interactions for comparison. At second order, the effective 
mass includes fc-mass and e-mass effects. 



gument k. As discussed in Sect. Ill Bt the antisym- 
metrized two-body interactions in the first- and second- 
order terms are given by vffl = (1 - Pi2)V\ ow k + 



V 3N /2 and V £ 



waves V t 



(2) 



(1 - Pl2) Mo 



(1) 



T^3N, with partial 

(2) 



Sll'J 



Viowk,swj + V 3 n,swj/4: and Vg U ',j 



Vlowk, Sll'J + ^3N, Sll'J /2. 

We solve the Dyson equation, Eq. (|28|) . self- 
consistently using the self-energies given by Eqs. ([50)1 
and (|31l) . In Fig. [51 we show the resulting effective mass 
at the Fermi surface, 



m*(fcp) 



k dk 



(32) 



k—k-e 



At the Hartree-Fock level, 3N contributions only change 
the effective mass marginally. Including second-order 
contributions leads to the typical enhancement of the ef- 
fective mass at the Fermi surface, and we find a larger 
impact of 3N forces for fcp > 1.3 fm - . 



C. Second order: energy per particle 

(2) (2) 

We include the second-order contributions E\ ' to E\' 
of Fig. 21 which are given by 

^ +3 N,e ff -J(nTr.J^)|(12|^)|34)| 2 
n kl n k2 (l - n k3 )(l - nkj , , 3 



£k 3 — £ k 4 



(2 7 r) d 5(k 1 +k 2 -k3-k 4 ). 



(33) 



As in the second-order self-energy, the antisymmetrized 

two-body interactions when evaluating contributions be- 

(2) 

yond the Hartree-Fock level are given by V£' = (1 - 
P12) Vi OVJ k + ^3N- The second-order calculations are car- 
ried out using the self-consistent single-particle energies 
determined by solving the Dyson equation, Eq. (|2"5)l . as 
discussed in Sect. MI 51 and the intermediate-state phase- 
space integrations are performed fully. Summin g oy er the 
spins and expanding in partial waves, we have |25| 

]T I (k SM S I K ( s 2) I k' SM' S ) 1 2 = J2 P l ( cos 0k,k<) 

S.Ms-.M's L 

x J2 Yl ^?i {l ~ l ' +T ~ v H-^) T+l ' +L c^c^ 



101'Q I'OIO 



x J {21 + 1)(2Z' + 1)(2Z + 1)(2Z' + 1) (2J + 1)(2J + 1) 



I S J 
J L V 



J S I' 
I L J 



(k\vg) u \k')(k'\V^\k) 



(2)_ 
Sl'U 1 



(1 - (-l) l+S+1 ) (1 - (-1) F + S+1 ) 



(34) 



where {. . .} denote 6j-symbols and Pl(cos9) are Legen- 
dre polynomials. Keeping only L — in Eq. (|34[) cor- 
responds to the angle-averaging approximation for the 
Pauli-blocking operator, but we keep all L ^ 6 for Viowfc 
and L ^ 4 for V 3 n- 

Our second-order results for the neutron matter en- 



ergy, i?NN+3N,cff = ^NN+3N,cff + ^NN+3N,cff ' are P re " 

sented in Fig. [3 The different contributions are listed 
in Table |U We observe that the cutoff dependence is 
reduced when going from first to second order. This is 
as expected based on the nuclear matter results 0, H| , 
but for neutron matter the cutoff dependence is signif- 
icantly weaker already at the Hartree-Fock level. The 
cutoff dependence increases with density and is less than 
1 MeV per particle for the densities studied in Fig. [7] 
over the cutoff range 1.8 far ^ A < 2.8 fm" 1 and 
2.0 fm -1 ^ A3NF ^ 2.5 fm -1 . This band sets the scale 
for omitted short-range many-body contributions and we 
discuss the theoretical uncertainties in the long-range 
parts in the following section. The weak cutoff depen- 
dence also demonstrates that the average momentum in 
the system (which is smaller than the Fermi momentum, 
because (pf) = 3/5 fcp) is well below the cutoff. 

Moreover, we have found that self-energy corrections 
to the neutron matter energy are practically negligible. 
The second-order energy with £ = is within 200 kcV 
of the self-consistent results shown in Fig. [7] In addi- 
tion, the second-order energy contributions are always 
below 1.3 MeV per particle in Table U except for the 
large cutoff A = 2.8 fm -1 cases. The second-order con- 
tributions practically only improve the cutoff indepen- 
dence of the results without changing the energy signifi- 
cantly. Moreover, for the lower cutoffs, the Hartree-Fock 
energies are already reliable. These findings combined 
suggest that neutron matter is perturbative at nuclear 



(i) 



E 



(2) 



9 




0.05 0.10 0.15 0.05 0.10 0.15 

P [fm" 3 ] p [fm" 3 ] 



FIG. 7: (Color online) Energy per particle E/N of neutron matter as a function of density p at the Hartree-Fock level (left) 
and including second-order contributions (right). The results are based on evolved N 3 LO NN potentials and N 2 LO 3N forces. 
Theoretical uncertainties are estimated by varying the NN cutoff (lines) and the 3N cutoff (band for fixed A = 2.0 fm -1 ). 
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-24.93 
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TABLE I: Contributions to the neutron matter energy due to the diagrams of Fig. H Results are given for Fermi momenta 
kF — 1.3, 1.5 and 1.7fm _1 and for different A/A3NF combinations. All energies are in MeV and /cf, A/A3NF are in fm -1 . 



densities. Therefore, we are confident that the P = 
approximation for Vsh is reliable, when evaluating the 
small second-order contributions, and that it is reason- 
able to neglect the residual 3N-3N diagram, E§ . 

D. Sensitivity to a uncertainties 

Next, we study the sensitivity of the second-order en- 
ergy to uncertainties in the coefficients that determine 
the long-range part of N 2 LO 3N forces. This provides an 
update for chiral potentials of the results of Ref. [25[ . The 



Ci coefficients relate 7rN, NN and 3N interactions, and 
the determination from 7rN scattering is, within errors, 
consistent with the extraction from NN waves. Present 
constraints for c\ and C3 are c\ — — 0.9^5 GeV -1 an d 
c 3 = -4.7tl;o GeV" 1 We note that at N 3 LO there 
are contributions that shift the Cj '10], and may lead to C3 
coefficients that are smaller in magnitude. In this study, 
we vary Cj only in 3N forces, because of lack of N 3 LO 
NN potentials that explore these Cj variations. However, 
based on the universality of V\ ovl k [1, E3 (starting from 
chiral potentials with two different c« sets [HI, [IB]), we 
do not expect large differences from varying c\ and C3 
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FIG. 8: (Color online) Theoretical uncertainties of the second- 
order energy with A/A3NP = 2.0 fm -1 as a function of density 
due to the uncertainties in the Ci and C3 coefficients of 3N 
forces. 

in NN interactions, where these variations are also ab- 
sorbed by higher-order contact interactions that have to 
be adjusted to reproduce NN scattering. 

In Fig. [8] we show that the theoretical uncertainties of 
the neutron matter energy are dominated by the uncer- 
tainties in the a coefficients, in particular the C3 part, 
compared to the uncertainties of the many-body calcula- 
tion or of neglected short-range many-body interactions 
probed by cutoff variations. The c\ and C3 variation leads 
to an energy uncertainty of ±1.5 MeV per particle at sat- 
uration density. Therefore, it is important to improve 
the constraints on c,. Figure [5] also shows that N 2 LO 3N 
forces provide a repulsive contribution to neutron matter 
at second order, even with this band. 

We compare our neutron matter energy with the Cj un- 
certainty band to other recent neutron matter results in 
Fig. H These include Monte-Carlo calculations (GFMC 
v 6 and vi HI and QMC S-wave @) and di-fermion EFT 
results [1] for lower densities, which are only sensitive to 
parts of NN interactions, at very low densities only to the 
neutron-neutron scattering length and effective range. In 
addition, we include in Fig. [5] the results of Akmal el 
al. [34[ based on the Argonne vis NN and Urbana IX 
3N potentials, where the repulsive 3N contributions arise 
from various terms. 



E. Symmetry energy 

The spread in the energy band of Fig. [5] corresponds 
to a range for the symmetry energy and its density de- 
pendence. In the following we explore this correlation 
and consequently also a possibility of using nuclear mat- 
ter properties to constrain some of the Cj couplings. As 



FIG. 9: (Color online) Comparison of the second-order energy 
with the d uncertainty band of Fig.[H]to other neutron matter 
results (see text for details). 



ci [GcV] 


c 3 [GeV] 


a 4 [MeV] 


po [MeVfm- 3 ] 


-0.81 


-3.2 


30.5 


2.0/2.0 


-0.81 


-5.7 


33.2 


2.8/2.8 


-0.7 


-3.2 


30.4 


2.0/2.0 


-1.4 


-5.7 


33.6 


2.8/2.9 



TABLE II: Symmetry energy £14 and the linear dependence 
of the symmetry energy po obtained from the energy band of 
Fig. [8] for different ci and C3 coefficients using ay = 16 MeV 
and K = 190/240 MeV. 



a function of asymmetry a = (p n — p P )/(p n + Pp), the 
energy of nuclear matter can be parameterized around 
saturation density as (see for example, Refs. [3ol [3rj|) 

^^ = -av + -^(p~p ) 2 + S 2 (p)a 2 + ... , (35) 

with the density-dependent symmetry energy £2 (p), 

S 2 (p) = a 4 + ^ (p - po) + . . . . (36) 
Po 

Here A is the nucleon number, ay the binding energy 
of symmetric nuclear matter at saturation density, K is 
the incompressibility, 04 the symmetry energy and po the 
linear dependence of the symmetry energy. 

Using the band of Fig. [8] over the density range 
0.13 fm -3 < p < 0.19 fm -3 with empirical values for the 
saturation point po — 0.16 fm~ 3 , ay = 16 MeV and in- 
compressibilities K = 190/240 MeV (which are in the 
range of empirical values [36( and were obtained for the 
same interactions in Ref . [8| ) , we find in Table |TT] that 
the symmetry energy ranges from 04 = (30.4 — 33.6) MeV 
and the linear dependence of the symmetry energy from 
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p Q = (2.0-2.9) MeVfm" 3 . The value for p is practically 
independent of the values of K. As shown in Table ITU the 
resulting range for the symmetry energy and its density 
dependence is set by the uncertainty of the dominant C3. 
Compared to the empirical range 04 = (25 — 35) MeV [3a |. 
the microscopic range of ~ 3 MeV is very useful and the 
comparison also suggests that smaller values in magni- 
tude for C3 may be somewhat favored. 



F. 1 So pairing gap 

Three-nucleon forces also impact superfluid pairing 
gaps or the anomalous self-energy. Here we will focus 
on the change of the neutron- neutron pairing gap A 
due to chiral 3N forces compared to the gap based on 
NN interactions. For simplicity, we neglect induced in- 
teractions and study the pairing gap at the BCS level. In 
this approximation, the contribution of 3N forces to the 
pairing interaction originates from two paired particles 
on the Fermi surface in back-to-back momentum states 
and the third particle is summed over occupied states in 
the superfluid BCS ground state. This means that the 
gap equation is given diagrammatically by 



(37) 



The kinematics in V3N is always P = for the interact- 
ing particle pair in the gap equation and therefore taking 
P = in T^3N is exact in this case. The only approx- 
imation using Vsn at the BCS level consists in having 
performed the sum over occupied states in the normal 
Fermi sea. However, the difference to summing over oc- 
cupied states in the BCS state is of higher order in A/ep 
(ef = kp/(2m) being the Fermi energy), which is small 
in neutron matter especially at the densities where 3N 
forces are effective. 

In the BCS approximation, the 1 So superfluid gap A(fc) 
is obtained by solving the gap equation, 



A(fc) 



1 



dpp 



V 1 f o (k, P )A(p) 
V(£ P -M) 2 + A 2 (p) 



(38) 



where the chemical potential /i = Ek F at the BCS level 
and the pairing interaction is given by the two-body 
interaction including V3N in the 1 Sq channel, so that 



(2) 



So -iowfc, 1 s + V Ws / 2 - The gap equation including 
3N forces, Eq. (|38|). can be derived explicitly by minimiz- 
ing the expectation value of the Hamiltonian in the BCS 
state H3- 

We consider two cases for the single-particle energies, 
a free spectrum — k 2 /(2m) and Hartree-Fock single- 
particle energies = 
the gap equation, Eq. 
suits for the neutron-neutron pairing gap on the Fermi 



k 2 /{2m) + £ (1) (fcV and solve 



as in Ref. 22]. Our re 
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FIG. 10: (Color online) The neutron-neutron So superfluid 
pairing gap on the Fermi surface A(kp) as a funtion of Fermi 
momentum &f. The results are based on evolved N 3 LO NN 
potentials and N 2 LO 3N forces with A/A3NF = 2.0 fm -1 using 
free and Hartree-Fock single-particle energies. 



surface A(fcp) are presented in Fig. [TO] as a function of 
the Fermi momentum. The inclusion of N 2 LO 3N forces 
leads to a reduction of the pairing gap for Fermi mo- 
menta fcp > 0.6 fm -1 , as expected since the channel 
of I^3N is repulsive in Fig. [3] This reduction agrees qual- 
itatively with results based on 3N potential models, see 
for example, Ref. [H)]. While the impact of 3N forces 
seems small in Fig. [101 for the densities, where the gap 
is decreasing, the reduction due to 3N forces is signifi- 
cant. Figure [10] also shows that the effects of 3N forces 
relative to NN interactions is very similar for free and 
Hartree-Fock single-particle energies. 



IV. CONCLUSIONS AND OUTLOOK 

In summary, we have shown that only the c\ and C3 
terms of the long-range 27r-exchange part of the leading 
chiral 3N forces contribute in neutron matter. Based on 
these parts and including all exchange terms, we derived 
density-dependent two-body interactions V3N for general 
momentum and spin configurations by summing the third 
particle over occupied states in the Fermi sea. The re- 
sulting V3N was found to be dominated by the repulsive 
central part. The comparison to the exact Hartree-Fock 
energy demonstrated that the P dependence is very weak 
and we have therefore taken P = in V3N- In addi- 
tion, the partial- wave matrix elements of our full V3N ap- 
proximately agree with the density-dependent two-body 
forces derived recently in the framework of in-medium 
chiral perturbation theory with certain approximations 
in Ref. M. 



12 



The density-dependent two-body interactions V3N cor- 
respond to the normal-ordered two-body part of 3N 
forces and we stressed that there are different normal 
ordering or symmetry factors when V3N is used at the 
Hartree-Fock, one- or two-body level. Moreover, it is im- 
portant that the summation of the third particle over 
occupied states is performed in the basis consistent with 
the one used in the many-body calculation. 

We have presented the first results for neutron matter 
based on chiral EFT interactions and including N 2 LO 
3N forces. The RG evolution of N 3 LO NN potentials to 
low momenta rendered the many-body calculation more 
controlled and our results for the energy suggest that 
neutron matter is perturbative at nuclear densities. This 
is based on small second-order contributions, self-energy 
corrections being negligible, and a generally weak cutoff 
dependence. We have found that N 2 LO 3N forces provide 
a repulsive contribution to the energy due to the repulsive 
central part in V3N. An important direction for future 
work is understanding the connection to the repulsive 
3N contributions to the oxygen isotopes discovered in 
Ref. [Ij. 

We have studied in detail the theoretical uncertainties 
of the neutron matter energy and found that the uncer- 
tainty in the C3 coefficient of 3N forces dominates com- 
pared to other many-body uncertainties. This resulted in 



an energy band in Fig. [9] Other recent neutron matter 
calculations were found to lie within this band. Moreover, 
the energy band provides microscopic constraints for the 
symmetry energy and its density dependence, and thus 
for the neutron skin in 208 Pb [42|. Finally, we have ob- 
tained a significant reduction of the 1 So superfluid pair- 
ing gap due to 3N forces for densities where the gap is 
decreasing. Our results show that chiral EFT and RG 
interactions can provide useful constraints for develop- 
ing a universal density functional based on microscopic 
interactions [3. IH l4l| . 
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